2.1 Artificial Neural Networks

Deep learning is an exciting invention that has risen in popularity in recent years, but it’s beginnings traced back to the 1950s when the earliest prototypes of artificial neural network algorithms were created. The algorithm is named so because it is inspired from our understanding at that time of how our biological brain responds to stimuli from sensory inputs. That is not to say that neural networks are valid representations of how our biological brain works - quite far from that! In fact, the over-sensationalization of neural network is in my opinion doing more harm to actual science than good.

To date, there are some school of thoughts that hypothesized about how our brain works. There’s the grandmother cell theory¹, which hypothesized that “neurons represent a complex but specific concept or object”. It was proposed by neuroscientists Jerry Lettvin, who suggested the notion of hyper-specific neurons that is activated by very specific prompts, such as upon hearing the mention of your grandmother, or seeing pictures of them. Advocates of this theory points to studies where participants were shown three to seven different “target” pictures of 132 photo subjects while having being observed - the researchers measured the electrical activity (activation) of the neurons connected to the electrodes and observe how neurons “almost ignored all other photos but fired steadily each time Aniston appeared on screen”².

Other theories assert that neurons are far simpler and that the representation of complex objects are distributed across many neurons. So while the human brain is extremely complex, with around 100 billion neurons and 1,000 trillion synaptic interconnections, each neuron rely on layers of electrochemical signaling and other biological machinery to process information.

Artificial neural networks are loosely derived from this second set of theories. A neural network architecture comprised of an arrangement of “layers” and “nodes” (representing neurons) such that “information” flows from one later and relayed to another.

Neuroscientists learn that³:
- Information from one neuron to another across a small gap called a synapse - One sending neuron can connect to several receiving neurons
- One receiving neuron can connect to several sending neurons

As one receiving neuron “receives” the signal sent by the sending neuron, the biochemical process help that neuron weight the signal according to its relative importance and as it turns out, this biological mechanism serves as a great template for an artificial neuron. The neurons are represented as nodes, and the strength between nodes are weighted to model the synapses in the biological brain.

As a motivational example, look at the following architecture:

By default, the neuralnet() function uses a backpropagation algorithm named rprop+. We are going to dive into the technical details of the backpropagation algorithm in later sections, but I want to take the opportunity to explain how this relates to “intelligence” as we understand it. A modern definition could be:

The ability to perceive information, and retain it as knowledge to be applied towards adaptive behaviors within an environment or context

So for a machine to fit that definition of intelligence, it has to not just “learn” the weights between these nodes, it has to learn how to adapt these weights - and that is achieved with what is termed as a cost function. The cost function determine how well the “learned weights” are by calculating the error given the paramaters. By means of analogy, you can imagine that a driver behind wheels having to approximate the optimal distance between his car and the car in front of him. Too much distance, and he may be driving too slow and causing unnecessary discomfort to other drivers around him. Too little distance and he may be accused of “tailgating” and put his own safety as well as others at risk. He / she has to learn the optimal distance through experience and “backpropagate” her visual feedback to correct the distance - when the car in front of her slows down, that visual feedback is passed and the “weights” are adjusted to accommodate that information. The intuition, or the sense of danger acts as a cost function, and she is trying to minimize the sense of danger or equivalently, the cost function.

Artificial intelligence imitates that objective when it is given a cost function because it now has a sense of badness and know how poorly the weights are in relation to the prediction error. A very simple cost function would be one iterative function (for-loop) where we iteratively make estimation, each time using a different combination of coefficients, and plot the error in each iteration by comparing our prediction to the ground truth.

Using the cost function, our artificial neural network can arrive at the optimal values for each coefficient that is being estimated and as a result of backpropagation, it updates all of the neuron paths (weights) in each iteration of feedback.

2.2 Neural Networks Architecture

Let’s plot our neural network again, and this time we’ll take the time to learn the topology of our network:

plot(nneu, rep = "best")

This is an architecturally simple neural network model, and while neural network architecture can vary in forms and design, each of them can be characterized by:
- Number of layers in the network
- Number of nodes in each layer
- Direction of signals in the network

When you wrap a neural network object, nn in plot(), the plot renders the network topology as well as the number of training steps that the neural network takes.

It also report the sum of squared errors, a measurement of error that all of you should be familiar by now. You can prevent the display of these error and steps by assigning an additional parameter parameter=FALSE in the plot() call.

Since the cost function we spoke of earlier is of such crucial importance to a neural network, let’s dive deeper into the mathematical details to strengthen and validate our understanding of this mechanism.

neuralnet(y ~ Var1 + Var2 + Var3, dat, hidden = 4, 
    rep = 5)

set.seed(1)
plot(neuralnet( y~Var1+Var2+Var3, dat, hidden=_______, rep=1), rep="best")

plot(nneu, rep = 4)

sse <- function(x) {
    return(sum((x - dat$y)^2))
}

length(nneu$net.result)

## [1] 5

### Write the code to inspect the prediction returned
### by any ONE of the 5 repetitions

# Write the code to compute your SSE for each of
# the 5 repetitions

set.seed(100)
nneu <- neuralnet(y ~ Var1 + Var2 + Var3, dat, hidden = 4, 
    rep = 5)
plot(nneu, rep = "best")

nneu_compt <- compute(nneu, dat[, 1:3], rep = 3)

data.frame(error_compute = nneu_compt$net.result, error_nn_3 = nneu$net.result[[3]])

##     error_compute      error_nn_3
## 1  0.007933439555  0.007933439555
## 2  0.985257670425  0.985257670425
## 3  1.002454293110  1.002454293110
## 4  0.002983716433  0.002983716433
## 5 -0.007984063057 -0.007984063057
## 6  1.003625540635  1.003625540635
## 7  0.992866239712  0.992866239712
## 8  0.010399328295  0.010399328295

1/2 * sum((dat$y - nneu_compt$net.result)^2)

## [1] 0.0002655641177

data.frame(computed = nneu$net.result[[3]], actual = dat$y)

##          computed actual
## 1  0.007933439555      0
## 2  0.985257670425      1
## 3  1.002454293110      1
## 4  0.002983716433      0
## 5 -0.007984063057      0
## 6  1.003625540635      1
## 7  0.992866239712      1
## 8  0.010399328295      0

plot(nneu, "best")

3.1 Backpropagation

So how does a neural network estimate the weights for each of the nodes? Well - remember that we already learned about a cost function which learns the error in each repetition. The key is to backpropagate these errors using an algorithm simply known as backpropagation.

Instead of approaching the estimation of weights using a more closed-form algebraic technique (as in OLS regression), the neural network model randomly initialize the weights and it learn the best “weights” by minimizing the cost function. For a regression problem, that is typically the sse that you’ve seen repeatedly up to this point of the workshop. For a classification problem, we can specify the cross-entropy error to be used (neuralnet(..., err.fct="ce")).

And recall for the error function that uses the sse: \(E = \frac{1}{2} \sum\limits^L_{l=1}\sum\limits^H_{h=1}(o_{lh}-y_{lh})^2\)

For an error function that uses the cross-entropy ce, it would be:
\(E = - \sum\limits^L_{l=1}\sum\limits^H_{h=1}(y_{lh}log(o_{lh}) + (1-y_{lh})log(1-o_{lh}))\)

Where:
- output o(x) is calculated for given inputs x and current weights
- \(l=1,...,L\) indexes the observations, i.e given input-output pairs
- \(h=1,...H\) are the output nodes

We’ve learned about the two functions in earlier workshops of the Machine Learning Specialization, but where the backpropagation starts to differ is that it modifies the weights of a neural network in order to find a local minimum of the error function with respect to the weights.

The learning of the error is the cost function of the neural network, and the feeding back of these error for the modification of weights is called backpropagation.

The de facto standard is through a general purpose optimization algorithm known as “gradient descent”.

library(animation)
## default example
ani.options(interval = 0.3, nmax = 50)
gd <- grad.desc()
# gd$par for the parameters gd$persp for
# perspective plot gd$persp(col = 'lightblue', phi
# = 30)

ani.options(interval = 0, nmax = 200)
f2 <- function(x, y) sin(1/2 * x^2 - 1/4 * y^2 + 3) * 
    cos(2 * x + 1 - exp(y))

gd2 <- grad.desc(f2, c(-2, -2, 2, 2), c(-1, 0.5), gamma = 0.1, 
    tol = 0.0001)

set.seed(100)
x <- as.matrix(dat[, 1:3])
x <- cbind(rep(1, nrow(x)), x)
y <- as.matrix(dat[, 4])
# initialize weights to 0.1
beta_gd <- matrix(0.01, nrow = ncol(x))
# learning rate
lrate <- 0.1

# repeat below for 4,000 iterations
for (i in 1:4000) {
    # cost function using vectorized operation
    cost <- (x %*% beta_gd) - y
    # calculate the gradient at that point
    delta <- (t(x) %*% cost) * (1/nrow(x))
    
    # move the guessed weights in opposite direction
    beta_gd <- beta_gd - (lrate * delta)
}
# linear regression w/ gradient descent
beta_gd

##               [,1]
##       5.1407142312
## Var1 -0.2432233805
## Var2  0.0684175269
## Var3  3.2402247446

# OLS closed-form solution
solve(t(x) %*% x) %*% t(x) %*% y

##               [,1]
##       5.1407142312
## Var1 -0.2432233805
## Var2  0.0684175269
## Var3  3.2402247445

summary(lm(y ~ ., dat))$coef[, 1]

##   (Intercept)          Var1          Var2 
##  5.1407142312 -0.2432233805  0.0684175269 
##          Var3 
##  3.2402247445

3.2 Neural Network for Classification

When we apply neural network as a classification algorithm, the network uses the logistic sigmoid function to squash the output of a linear regression to a value within the bounds of 0 to 1. This, you may recall, is what you’ve learned in the Classification 1 workshop. Allow me to help jog your memory:

• A logistic regression is a linear regression technique modified to solve binary classification problems
  
• The “activation” function of a classification-purposed neural network is the logistic function (transform the output of a linear regression to a value within the bounds of 0 and 1)

I will spare you the detailed explanation of the logistic regression but the following code is identical to the one you’ll find in your Classification 1 coursebook. Exponenting something would always make it a positive value, and then to make the final value lesser than 1, we divide it by a number greater than itself (+1 arbitrarily):

sigmoid <- function(x) {
    exp(x)/(1 + exp(x))
}
# equivalent: sigmoid <- function(x){ 1/(1+exp(-x))
# }

sigmoid(-5)

## [1] 0.006692850924

Let’s confirm that the function we wrote above does resemble a sigmoid curve:

curve(sigmoid, xlim = c(-10, 10))

Let’s also create our cross-entropy loss function:

ce_manual <- function(x, y) {
    -sum((y * log(x) + (1 - y) * log(1 - x)), na.rm = T)
}

This of course, is just the code implementation of the ce cross-entropy error that we learned in earlier chapters: \(E = - \sum\limits^L_{l=1}\sum\limits^H_{h=1}(y_{lh}log(o_{lh}) + (1-y_{lh})log(1-o_{lh}))\)

We’ll generate some data, this time with y being a binary class:

set.seed(100)
y <- c(0, 1, 1, 0)
dat <- data.frame(expand.grid(c(0, 1), c(0, 1), c(1, 
    2)), y)
head(dat)

##   Var1 Var2 Var3 y
## 1    0    0    1 0
## 2    1    0    1 1
## 3    0    1    1 1
## 4    1    1    1 0
## 5    0    0    2 0
## 6    1    0    2 1

And we can write our neural network algorithm by combining the sigmoid() function we wrote earlier with the backpropagation technique. We will also compute the cross-entropy error to see if our manually-written algorithm can rival the neuralnet() algorithm in performance:

neuralnet_logit <- function(x, y, epochs, lrate = 0.005) {
    x <- cbind(1, x)
    beta_hat <- matrix(0.1, nrow = ncol(x))
    for (i in 1:epochs) {
        error <- sigmoid(x %*% beta_hat) - y
        # update weights with gradient descent
        delta <- t(x) %*% as.matrix(error, ncol = nrow(x)) * 
            (1/nrow(x))
        beta_hat <- beta_hat - (lrate * delta)
    }
    print(ce_manual(sigmoid(x %*% beta_hat), y))
}

And with that we’re ready to apply the algorithm we just wrote to the data, with 500 iterations of the update rule. Here I’m introducing a new terminology called epochs, which is a common term in neural network terminology. One epoch consist of one round of forward pass and one round of backward propagation of the error. Another way to put it:

An epoch is a measure of the number of times all of the training vectors are used once to update the weights. For batch training all of the training samples pass through the learning algorithm simultaneously in one epoch before weights are updated.

neuralnet_logit(as.matrix(dat[, 1:3]), dat$y, epochs = 500)

## [1] 5.54745663

epoch is not to be confused with the rep parameter, with the latter actually referring to the number of repetitions for the neural network’s training. The log-likelihood (or cross-entropy) error is 5.547. As a reminder, a perfect binary classifier would have an error of 0:

ce_manual(dat$y, dat$y)

## [1] 0

Let’s now use the neuralnet package: the parameters are largely similar, except we override the err.fct with ce and set linear.output to FALSE since we are configuring the neural network to use a logistic function instead of the linear one earlier:

set.seed(100)
nn_logit <- neuralnet(y ~ Var1 + Var2 + Var3, dat, 
    hidden = 0, err.fct = "ce", linear.output = F, 
    rep = 500)
ce_manual(nn_logit$net.result[[1]], dat$y)

## [1] 5.545689639

Discussion: The local minima that the function arrive at doesn’t seem to be a lot better than ours. How would you improve the neural network model? - [ ] Set higher amount of rep for the neural network - [ ] Change the seed to have a different set of initial weights
- [ ] Change err.fct from ce to sse
- [ ] Add a hidden layer to the neural network in the hope that the increased complexity would increase its ability to model the pattern within the data

If you pick the fourth option and plot the best repetition, you’ll arrive at the following neural network plot:

set.seed(100)
nn_logit <- neuralnet(y ~ Var1 + Var2 + Var3, dat, 
    hidden = 3, err.fct = "ce", linear.output = F, 
    rep = 500)
plot(nn_logit, rep = "best")

We can manually confirmed that the error is in fact the ce error that we define by hand earlier. The function we wrote is ce_manual() so let’s verify that:

nn_logit_rs <- compute(nn_logit, dat[, 1:3], rep = which.min(nn_logit$result.matrix[1, 
    ]))$net.result
ce_manual(nn_logit_rs, dat$y)

## [1] 0.01875205382

While a multilayer perceptron may seem a little complicated (and that’s because it is complicated), the interpretation and mathematics are not any more difficult. Whether for regression or for classification, the difference is in the activation function (one is linear, while the other is sigmoid) and the cost function (one is sse and the other is ce). Consider the following diagram¹¹:

To compute the value of our output unit \(a_5\) using the multi-layer perceptron:

\(a_5 = g(w_{3,5} \cdot a_3 + w_{4,5} \cdot a_4)\)
\(a_5 = g(w_{3,5} \cdot g(w_{1,3} \cdot a_1 + w_{2,3} \cdot a_2) + w_{4,5} \cdot g(w_{1,4} \cdot a_1 + w_{2,4} \cdot a_2) )\)

par(mfrow = c(2, 2))
curve(plogis, from = -5, to = 5, main = "Logistic Regression")
curve(tanh(x), xlim = c(-5, 5), main = "Hyperbolic Tangent")
curve(exp(x)/sum(exp(x)), xlim = c(-5, 5), main = "Softmax Function")
curve(ifelse(x >= 0, x, 0), xlim = c(-5, 5), main = "Rectified Linear Unit")

6.1 [Optional] Neural Interpretation Diagram (NID)

Olden & Jackson addresses these issues concerning neural network visualization and put forth the Neural Interpretation Diagram¹⁶, which is then implemented in the NeuralNetTools library through the plotnet() function.

In the following, we plotted the same network we have created in earlier exercises (nneu) but as a neural interpretation diagram as in Ozesmi and Ozesmi (1999)¹⁷. By default, this diagram color positive weights between neurons as black lines and negative weights as gray lines:

suppressPackageStartupMessages(library(NeuralNetTools))
plotnet(nneu)

We can additionally specify our preferred colors for the circle (circle_col), border (bord_col), line colors for positive and negative weights (pos_col and neg_col respectively), and even adjust the color transparency:

plotnet(nneu, circle_col = "dodgerblue1", bord_col = "black", 
    pos_col = "darkblue", neg_col = "black", alpha_val = 0.6)

Applying the above set of parameters in our plotnet() function call but on the bank direct marketing neural network model:

plotnet(bankmar_nn2, circle_col = "dodgerblue1", bord_col = "dodgerblue1", 
    pos_col = "darkblue", neg_col = "black", alpha_val = 0.6)

7.1 Thinking about Parallelism

MXNet supports training with multiple CPUs and GPUs, which may be located on different physical machines. It is thus well-suited for learning tasks that involve a large amount of data and when the procedure can be parallelized. As it turns out, a deep learning network is one such model that allows the computation of parameters to be handled by different devices. A machine learning algorithm that is Bayesian in nature is not as suited from parallelism. In contrast to model parallelism is the notion of data parallelism, where each device is assigned an equal proportion of the data and the output from each device is aggregated to produce the final model.

By default, MXNet partitions a data batch evenly among the available GPUs: assume a batch size b and assume there are k GPUs, then in one iteration each GPU will perform the forward and backward on b/k examples. The gradients are then summed over all GPUs before updating the model.

And on model parallelism, its documentation says:

Model parallelism in deep learning was first proposed for the extraordinarily large convolutional layer in GoogleNet. From this implementation, we take the idea of placing each layer on a separate GPU. Using model parallelism in such a layer-wise fashion provides the benefit that no GPU has to maintain all of the model parameters in memory.

The official documentation provides specific installation instructions. Alternatively, you can refer to the README on MXNet’s GitHub page.

# cran <- getOption('repos') cran['dmlc'] <-
# 'https://apache-mxnet.s3-accelerate.dualstack.amazonaws.com/R/CRAN/'
# options(repos = cran) install.packages('mxnet')
# if on R 3.5 and the above doesnt work, then try:
# install.packages('https://s3.ca-central-1.amazonaws.com/jeremiedb/share/mxnet/CPU/mxnet.zip',
# repos = NULL)
library(mxnet)

Now that you have MXNet on your operating system, we’re ready for digit classification task. A classic introduction to deep learning or computer vision is the MNIST dataset, which stands for “Modified National Institute of Standards and Technology”. The goal is to idenfity digits from a dataset of tens of thousands of handwritten images, and the dataset is attributed to Yann LeCun¹⁸.

We’ll read the data into our environment:

mnist_train <- read.csv("data_input/mnist/train.csv")
mnist_test <- read.csv("data_input/mnist/test.csv")

Tip: If the loading is too slow on your machine, consider using library(readr) and then replace read.csv with read_csv.

Discussion:
1. What do you think the original image dimension is? (Tip: all images are squares)

1. 72px x 72px
   
2. 36px x 36px
   
3. 28px x 28px

2. The target variable is label, are the distribution of our class labels equally distributed? With image processing you want a fairly equal distribution so each class is sufficiently exposed to the neural network.

Our first column is label, and the rest are values between 0 and 255 representing grayscale numbers of handwritten digits in the corresponding pixel.

colnames(mnist_train)[c(1:5, 780:784)]

##  [1] "label"    "pixel0"   "pixel1"   "pixel2"  
##  [5] "pixel3"   "pixel778" "pixel779" "pixel780"
##  [9] "pixel781" "pixel782"

Sure enough, the distribution of our class labels do seem to be spread out quite evenly. While you may be tempted to use pie() to represent the distribution, you will want to resist that temptation - pie charts are generally frowned upon in the scientific community because the human eyes are not good at discerning the distribution from a pie chart.

tr_lab <- table(mnist_train$label)
barplot(tr_lab, main = "Distribution of Digits in Training Sample", 
    col = gray.colors(10))

As it was the case earlier, our neural network implementation rely on matrix algebra and we’ll want to represent our data as matrix objects. Be careful that you do not include the first column in your matrix creation (use: 2:ncol() or -1 to exclude the first column). The matrix would be a 28 x 28 and should be numeric (we’ll use as.numeric as an additional assurance).

Finally, because the matrix to be passed when creating our image later will be “flipped” vertically, we shall flip it back and fortunately R has a built-in way of handling that: the rev() function returned a reversed version of its argument, such that:

mat <- matrix(c(2:10), ncol = 3)
mat

##      [,1] [,2] [,3]
## [1,]    2    5    8
## [2,]    3    6    9
## [3,]    4    7   10

apply(mat, 2, rev)

##      [,1] [,2] [,3]
## [1,]    4    7   10
## [2,]    3    6    9
## [3,]    2    5    8

Putting it together, we now convert the 784 predictor (columns) of our first row into a matrix of 28 x 28 shape, ensure that the elements are numeric (doesn’t make sense to perform matrix computation on non-numeric values) and finally reverse the order of the rows in our matrix:

m1 <- matrix(mnist_train[1, 2:ncol(mnist_train)], nrow = 28, 
    ncol = 28)
m1 <- apply(m1, 2, as.numeric)
m1 <- apply(m1, 2, rev)

We should expect m1 to resemble a numeric “1” when we plot it using image() later.

head(mnist_train[1, 1])

## [1] 1

We’ve learned about the technique to display a grid of colored “pixels” with the colors corresponding to the values in a matrix in the Unsupervised Machine Learning course under the “eigenfaces” chapter. Here we’re using the same technique to render our m1 matrix as an image grid. By default, image() uses col.heat.colors(12) but that’s not something we’re going to worry about right now.

We added a text onto that image in the second line of code, specifying mnist_train[1,1] so the “ground truth” for row 1 of the original data would be displayed too as a reference:

image(1:28, 1:28, m1)
text(2, 2, col = "white", cex = 1.2, mnist_train[1, 
    1])

With that sanity test, we can now wrap a for-loop function over the same code so that each row of our data (for i in 1:nrow(data)) is visualized using the same image() function. This time, within the image() function call, we’ll specify grey.colors(255) instead of the default heat.colors(12). The decision to use 255 shades of gray isn’t exactly arbitrary - it is because in the original data our values actually assume a range of 0 to 255 as well:

range(mnist_train[1, 2:ncol(mnist_train)])

## [1]   0 255

Another detail is the use of sqrt(ncol(x)) to set the dimensions of our resulting image. We know our image is square (same number of row and columns) so this trick works. Instead of hard-coding a text value at the corner of each image, we will use i as the placeholder as well (mnist_train[i,1]). This two techniques combined help make sure our code will run is as general an environment as possible (no hard-coding is always better), where any data frame can be accepted and the resulting images can be plotted:

vizTrain <- function(input) {
    
    dimmax <- sqrt(ncol(mnist_train[, -1]))
    
    dimn <- ceiling(sqrt(nrow(input)))
    par(mfrow = c(dimn, dimn), mar = c(0.1, 0.1, 0.1, 
        0.1))
    
    for (i in 1:nrow(input)) {
        m1 <- matrix(input[i, 2:ncol(input)], nrow = dimmax, 
            byrow = T)
        m1 <- apply(m1, 2, as.numeric)
        m1 <- t(apply(m1, 2, rev))
        
        image(1:dimmax, 1:dimmax, m1, col = grey.colors(255), 
            xaxt = "n", yaxt = "n")
        text(2, 20, col = "white", cex = 1.2, mnist_train[i, 
            1])
    }
    
}

Let’s do another sanity test to make sure the function we created, vizTrain works as expected. Can you modify the code below to plot the first 36 digits?

vizTrain(mnist_train[1:25, ])

mxnet has its own “grammar” in defining a neural network’s architecture. Specifically, it made use of what it termed as the symbolic API, starting with a placeholder data symbol.

In the following code, we’re assembling a network with the following architecture: data --> 128 units ReLU --> 64 units ReLU --> 10 units output softMax

# notice how our layers are fully-connected and
# cascading
m1.data <- mx.symbol.Variable("data")

m1.fc1 <- mx.symbol.FullyConnected(m1.data, name = "fc1", 
    num_hidden = 128)
m1.act1 <- mx.symbol.Activation(m1.fc1, name = "activation1", 
    act_type = "relu")

m1.fc2 <- mx.symbol.FullyConnected(m1.act1, name = "fc2", 
    num_hidden = 64)
m1.act2 <- mx.symbol.Activation(m1.fc2, name = "activation2", 
    act_type = "relu")

m1.fc3 <- mx.symbol.FullyConnected(m1.act2, name = "fc3", 
    num_hidden = 10)
m1.softmax <- mx.symbol.SoftmaxOutput(m1.fc3, name = "softMax")

Many of the parameters are up to tuning, but what you absolutely want to do is to make sure that your last layers have 10 units, corresponding to the 10 output class.

As a summary:
- mx.symbol.Variable is where we define variables
- mx.symbol.FullyConnected is where set configure the input layer
- mx.symbol.Activation is where we set our activation function (sigmoid, relu, tanh or softmax typically)
- mx.symbol.Output is where we configure the final output layer

The reason why I used a softMax in the final layer is something I discussed in Remarks 2 in the earlier chapter: it is convenient for a multi-class MLP network (which is what a digit recognition is) because the sum of probabilities of all possible class will be 1. If this is hard to intuit about, just recall the formula to plot the softMax curve:

curve(exp(x)/sum(exp(x)), xlim = c(-5, 5), main = "Softmax Function")

The denominator is simply the sum, so the sum of exp(x) will result in 1.

c(4, 8, 1)/sum(c(4, 8, 1))

## [1] 0.30769230769 0.61538461538 0.07692307692

So we have a neural network created using the mxnet specification model, data --> 128 units ReLU --> 64 units ReLU --> 10 units output softMax

But that’s not very friendly on the eyes; MXNet uses the graph.viz to graph the architecture of our network and allow us to visualize an otherwise complex (and messy) network:

graph.viz(m1.softmax)

Now that we have our data and our network ready, the rest is something I believe all of you should be familiar with. We’ll go ahead and prepare the train-test split just as we’ve always done. First we’ll convert the data frames back into matrix:

train <- data.matrix(mnist_train)
test <- data.matrix(mnist_test)

The training data would

# remove 'label' train_x: result in 42000 x 784
# matrix range 0 to 255
train_x <- train[, -1]
train_y <- train[, 1]

# train_x: result in 784 x 42000 matrix range 0 to
# 1 (divide by 255) each of the 784 is treated as
# input (predictors), hence the transpose
train_x <- t(train_x/255)
test_x <- t(test/255)

To create a MXNet feedforward neural network, we’ll pass in the symbolic network specification we created above, our X (input) and y (output). Most of the other variables are optional and have sensible default values. num.round, for example by default uses 10: that’s how we specify the number of iterations over training data to train the model. I’ve increased this value to 40;

array.layout is set to auto by default. For a matrix we use rowmajor when our dimensions are c(example, features) and colmajor when they are c(features, examples). When it’s auto, then it will try to auto-detect the layout by matching the feature size. We also specified “accuracy” as our performance metric, and finally we specified a callback to log the time our neural network uses:

log <- mx.metric.logger$new()
startime <- proc.time() 
mx.set.seed(0)

m1 <- mx.model.FeedForward.create(m1.softmax,  #the network configuration made above
                                     X = train_x, #input (predictors),
                                     y = train_y, #the labels
                                     ctx = mx.cpu(),
                                     num.round = 40,  # 10 get us ~0.93 accuracy
                                     array.batch.size = 80,
                                     momentum = 0.95,
                                     array.layout="colmajor",
                                     learning.rate = 0.001,
                                     eval.metric = mx.metric.accuracy,
                                     epoch.end.callback = mx.callback.log.train.metric(1,log)
)

## Start training with 1 devices

## [1] Train-accuracy=0.110642857704134

## [2] Train-accuracy=0.111523810018386

## [3] Train-accuracy=0.111523810018386

## [4] Train-accuracy=0.111523810018386

## [5] Train-accuracy=0.156023809366993

## [6] Train-accuracy=0.24204761810246

## [7] Train-accuracy=0.533166666541781

## [8] Train-accuracy=0.761976190748669

## [9] Train-accuracy=0.81309523945763

## [10] Train-accuracy=0.835142857233683

## [11] Train-accuracy=0.850523810273125

## [12] Train-accuracy=0.862476190044766

## [13] Train-accuracy=0.876404762381599

## [14] Train-accuracy=0.889880952608018

## [15] Train-accuracy=0.901357142925262

## [16] Train-accuracy=0.909285715307508

## [17] Train-accuracy=0.915857142266773

## [18] Train-accuracy=0.921499999704815

## [19] Train-accuracy=0.92645238115674

## [20] Train-accuracy=0.9316428565979

## [21] Train-accuracy=0.936476189749581

## [22] Train-accuracy=0.940380950655256

## [23] Train-accuracy=0.944238094375247

## [24] Train-accuracy=0.947380950905028

## [25] Train-accuracy=0.950261904171535

## [26] Train-accuracy=0.95276190394447

## [27] Train-accuracy=0.954928570815495

## [28] Train-accuracy=0.95697618995394

## [29] Train-accuracy=0.958619046665373

## [30] Train-accuracy=0.960452380066826

## [31] Train-accuracy=0.961952380112239

## [32] Train-accuracy=0.9635714277767

## [33] Train-accuracy=0.965380952812376

## [34] Train-accuracy=0.967071429093679

## [35] Train-accuracy=0.968142857892173

## [36] Train-accuracy=0.969523810659136

## [37] Train-accuracy=0.970619048731668

## [38] Train-accuracy=0.971571430478777

## [39] Train-accuracy=0.972333335308802

## [40] Train-accuracy=0.973571431069147

print(paste("Training took:", round((proc.time() - startime)[3],2),"seconds"))

## [1] "Training took: 173.87 seconds"

In less than 100 seconds we’ve trained a neural network that classify handwritten digits with 97.4% accuracy! In the network configuration, we specified momentum, a concept strongly related to the gradient descent algorithm we spoke about earlier.

Recall from earlier lessons that gradient descent is a greedy algorithm that at each step selects the steepest descent at the current point as the direction of advancement. This means that the algorithm may in many cases, instead of descending towards the center directly would zig-zag around because of the noise in our sample size when calculating the gradient:

In the case above, the curvature in the direction of the x2 axis is far larger than that of the x1, causing the gradient descent to advance more aggressively towards x2. The momentum algorithm addresses that problem by introducing a parameter \(\gamma\) ranging from 0 to 1 to the descent, and is so doing increases the size of steps taken towards the minimum (jumping out from a local minima).

Say our first update to the weights is \(\theta_1\) and the second update without momentum is \(\theta_2\). However, with momentum, the second update would instead have been \(\theta_2 + \gamma \cdot \theta_1\). For yet another update, it would have been \(\theta_3 + \gamma \cdot \theta_2 + \gamma^2 \cdot \theta_1\) and so on. The objective is to avoid the algorithm getting stuck in a local minimum when the error function isn’t a smooth but consist of several local minima:

By introducing momentum into the gradient descent, it also pushes the algorithm to converge faster to the global minimum:
Illustration: Momentum from Sratch¹⁹

By default, mx.model.FeedForward.create assumes verbose=TRUE and prints out the information of each iteration during training. We can suppress that to false, and then plot the accuracy at each iteration using plot(log$train) instead. This is useful in production environment when we’re typically running 500 iterations or more and do not require information being printed to the console. For this technique to work, you want to initialize the mxnet logger by calling mx.metric.logger$new() before training your deep learning network.

plot(log$train, type = "l", xlab = "Iteration", ylab = "Accuracy")

7.2 Making Predictions

Up to this point we haven’t tested our model on the test dataset. The process here is similar to what you’ve learned in the past: using predict() and specify that the data is in the colmajor format (c(features, examples)). If you want to confirm these, dim(test_x) is a quick and easy test.

The predict() function will return a prediction in the same colmajor format, so we’ll expect 10 rows and 28000 columns, each value corresponding to the probability of an example being in a particular class. If you rather see each example as a row, and probability of it being a digit as a column ,then you can use t() to transpose the result:

m1_preds <- predict(m1, test_x, array.layout = "colmajor")
t(round(m1_preds[, 1:5], 2))

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.00    0 1.00 0.00 0.00 0.00 0.00 0.00 0.00
## [2,] 1.00    0 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## [3,] 0.00    0 0.00 0.00 0.10 0.00 0.00 0.00 0.04
## [4,] 0.07    0 0.02 0.00 0.49 0.01 0.01 0.04 0.04
## [5,] 0.00    0 0.01 0.99 0.00 0.00 0.00 0.00 0.00
##      [,10]
## [1,]  0.00
## [2,]  0.00
## [3,]  0.84
## [4,]  0.33
## [5,]  0.00

Notice that our neural network was very sure in its prediction for the first three examples (2, 0, 9 respectively), but less-certain on the forth one even though it predicted a 9.

Let’s transpose that back to it’s output format and print the first 10 predictions by taking the column which has the largest probability using max.col(). We also added -1 because the index should correspond to a digit smaller than itself (so 4th column correspond to digit 3, and 10th column correspond to digit 9):

m1_preds_result <- max.col(t(m1_preds)) - 1
m1_preds_result[1:10]

##  [1] 2 0 9 4 3 7 0 3 0 3

We can also modify the plotting function we’ve written above, to plot the prediction from m1_preds_result right on each training example. This is roughly the same code, with the difference being the text(..., preds[i]) which should print out the prediction:

plotResults <- function(images, preds) {
    
    x <- ceiling(sqrt(length(images)))
    par(mfrow = c(x, x), mar = c(0.1, 0.1, 0.1, 0.1))
    
    for (i in images) {
        m <- matrix(test[i, ], nrow = 28, byrow = TRUE)
        m <- apply(m, 2, rev)
        image(t(m), col = grey.colors(255), axes = FALSE)
        text(0.05, 0.1, col = "red", cex = 1.2, preds[i])
    }
    
}

Here let’s put our visualization function to work. Pay special attention to the digit that our neural network wasn’t sure about (the fourth digit it says was a 9 but was only 64% sure whereas 95%-100% certain on pretty much every other prediction):

plotResults(1:49, m1_preds_result)